library(knitr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(raster)
library(rasterVis)
library(scales)
getData("ISO3")%>%
as.data.frame%>%
filter(NAME=="Bangladesh")
## ISO3 NAME
## 1 BGD Bangladesh
Often good idea to keep data in separate folder. You will need to edit this for your machine!
datadir="~/GoogleDrive/Work/courses/2015_UB503/SpatialR/data"
Download country border.
bgd=getData('GADM', country='BGD', level=0,path = datadir)
plot(bgd)
SRTM Elevation data with getData() as 5deg tiles.
bgdc=gCentroid(bgd)%>%coordinates()
dem1=getData("SRTM",lat=bgdc[2],lon=bgdc[1],path=datadir)
plot(dem1)
plot(bgd,add=T)
Download the remaining necessary tiles
dem2=getData("SRTM",lat=23.7,lon=85,path=datadir)
Use merge() to join two aligned rasters (origin, resolution, and projection). Or mosaic() combines with a function.
dem=merge(dem1,dem2)
plot(dem)
plot(bgd,add=T)
Beware of massive temporary files!
inMemory(dem)
## [1] FALSE
dem@file@name
## [1] "/private/var/folders/lg/gz_jcfk5617dlpzdtg23x3rh0000gn/T/R_raster_adamw/r_tmp_2015-10-19_114157_62728_62244.grd"
file.size(sub("grd","gri",dem@file@name))*1e-6
## [1] 144
showTmpFiles()
## r_tmp_2015-10-19_105338_61539_09013
## r_tmp_2015-10-19_105340_61539_06688
## r_tmp_2015-10-19_105809_61709_12650
## r_tmp_2015-10-19_105934_61709_54451
## r_tmp_2015-10-19_105939_61709_84199
## r_tmp_2015-10-19_110857_61940_90130
## r_tmp_2015-10-19_110900_61940_66881
## r_tmp_2015-10-19_111013_61940_26505
## r_tmp_2015-10-19_111142_61940_44518
## r_tmp_2015-10-19_111147_61940_41991
## r_tmp_2015-10-19_111532_54312_02140
## r_tmp_2015-10-19_113231_62659_02529
## r_tmp_2015-10-19_113242_62673_02529
## r_tmp_2015-10-19_113257_62686_02529
## r_tmp_2015-10-19_113326_62711_02529
## r_tmp_2015-10-19_114157_62728_62244
rasterOptions()
## format : raster
## datatype : FLT8S
## overwrite : FALSE
## progress : none
## timer : FALSE
## chunksize : 1e+07
## maxmemory : 1e+08
## tmpdir : /var/folders/lg/gz_jcfk5617dlpzdtg23x3rh0000gn/T/R_raster_adamw/
## tmptime : 168
## setfileext : TRUE
## tolerance : 0.1
## standardnames : TRUE
## warn depracat.: TRUE
## header : none
Set with rasterOptions(tmpdir = "/tmp")
Saving raster to file: two options
Save while creating
dem=merge(dem1,dem2,filename=file.path(datadir,"dem.tif"))
Or after
writeRaster(dem, filename = file.path(datadir,"dem.tif"))
| Filetype | Long name | Default extension | Multiband support |
|---|---|---|---|
| raster | ‘Native’ raster package format | .grd | Yes |
| ascii | ESRI Ascii | .asc | No |
| SAGA | SAGA GIS | .sdat | No |
| IDRISI | IDRISI | .rst | No |
| CDF | netCDF (requires ncdf) |
.nc | Yes |
| GTiff | GeoTiff (requires rgdal) | .tif | Yes |
| ENVI | ENVI .hdr Labelled | .envi | Yes |
| EHdr | ESRI .hdr Labelled | .bil | Yes |
| HFA | Erdas Imagine Images (.img) | .img | Yes |
rgdal does even more…
dem=crop(dem,bgd,filename=file.path(datadir,"dem_bgd.tif"),overwrite=T)
plot(dem); plot(bgd,add=T)
gplot(dem,max=1e5)+geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("red","yellow","grey30","grey20","grey10"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e3)))+
coord_equal(ylim=c(-21,25))+
geom_path(data=fortify(bgd),
aes(x=long,y=lat,order=order,group=group),size=.5)
## Regions defined for each Polygons
terrain() options:
Use a smaller region:
reg1=crop(dem1,extent(93.8,94,21.05,21.15))
plot(reg1)
The terrain indices are according to Wilson et al. (2007), as in gdaldem.
slope=terrain(reg1,opt="slope",unit="degrees")
plot(slope)
aspect=terrain(reg1,opt="aspect",unit="degrees")
plot(aspect)
Difference between the value of a cell and the mean value of its 8 surrounding cells.
tpi=terrain(reg1,opt="TPI")
gplot(tpi,max=1e6)+geom_tile(aes(fill=value))+
scale_fill_gradient2(low="blue",high="red",midpoint=0)+
coord_equal()
Negative values indicate valleys, near zero flat or mid-slope, and positive ridge and hill tops
plot() to:
Hint: use transparent to plot a transparent pixel.
Extract peaks/ridges and valleys:
plot(reg1)
plot(tpi>15,col=c("transparent","red"),add=T,legend=F)
plot(tpi<(-12),col=c("transparent","blue"),add=T,legend=F)
Mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.
tri=terrain(reg1,opt="TRI")
plot(tri)
Difference between the maximum and the minimum value of a cell and its 8 surrounding cells.
rough=terrain(reg1,opt="roughness")
plot(rough)
Compute from slope and aspect (in radians). Often used as a backdrop for another semi-transparent layer.
hs=hillShade(slope*pi/180,aspect*pi/180)
plot(hs, col=grey(0:100/100), legend=FALSE)
plot(reg1, col=terrain.colors(25, alpha=0.5), add=TRUE)
Flow direction (of water), i.e. the direction of the greatest drop in elevation (or the smallest rise if all neighbors are higher).
Encoded as powers of 2 (0 to 7). The cell to the right of the focal cell ‘x’ is 1, the one below that is 2, and so on:
| 32 | 64 | 128 |
|---|---|---|
| 16 | x | 1 |
| 8 | 4 | 2 |
flowdir=terrain(reg1,opt="flowdir")
plot(flowdir)
Much more powerful hydrologic modeling in GRASS GIS
slr=data.frame(year=2100,
scenario=c("RCP2.6","RCP4.5","RCP6.0","RCP8.5"),
low=c(0.26,0.32,0.33,0.53),
high=c(0.54,0.62,0.62,0.97))
kable(slr)
| year | scenario | low | high |
|---|---|---|---|
| 2100 | RCP2.6 | 0.26 | 0.54 |
| 2100 | RCP4.5 | 0.32 | 0.62 |
| 2100 | RCP6.0 | 0.33 | 0.62 |
| 2100 | RCP8.5 | 0.53 | 0.97 |
Range from 2.5-10m in Bangladesh since 1960 Karim & Mimura, 2008.
ss=c(2.5,10)
1st Question: How much area is likely to be flooded by rising sea levels?
WGS84 data is unprojected, must account for cell area (in km^2)…
area=area(dem)
plot(area)
Steps:
cellStats() to calculate potentially flooded areas.flood1=dem<=2.76
flood2=dem<=10.97
plot(flood2,col=c("transparent","darkred"))
plot(flood1,col=c("transparent","red"),add=T)
flood1a=flood1*area
flood2a=flood2*area
cellStats(flood1a,sum)
## [1] 6358.114
cellStats(flood2a,sum)
## [1] 68051.34
Socioeconomic Data and Applications Center (SEDAC) http://sedac.ciesin.columbia.edu
Data not available for direct download (e.g. download.file())
A nicer plot…
gplot(pop,max=1e6)+geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("grey90","grey60","darkblue","blue","red"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e5)))+
coord_equal()
dem objectpop2=pop%>%
crop(dem)
gplot(pop2,max=1e6)+geom_tile(aes(fill=value))+
scale_fill_gradientn(colours=c("grey90","grey60","darkblue","blue","red"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e5)))+
coord_equal()
Assume equal density within each grid cell and resample
pop3=pop2%>%
resample(dem,method="bilinear")
gplot(pop3,max=1e6)+geom_tile(aes(fill=value))+
scale_fill_gradientn(colours=c("grey90","grey60","darkblue","blue","red"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e5)))+
coord_equal()
Or resample elevation to resolution of population:
res(pop2)/res(dem)
demc=dem%>%
aggregate(fact=50,fun=min,expand=T)%>%
resample(pop2,method="bilinear")
How many people are likely to be displaced?
Steps:
flood2) x population density x areacellStats()flood2floodpop2=flood2a*pop3
cellStats(floodpop2,sum)
## [1] 82577156
Number of potentially affected people across the region.
gplot(floodpop2,max=1e6)+geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("grey90","grey60","darkblue","blue","red"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e4)))+
coord_equal()
distance() calculates distances for all cells that are NA to the nearest cell that is not NA.
popcenter=pop2>5000
popcenter=mask(popcenter,popcenter,maskvalue=0)
plot(popcenter,col="red",legend=F)
In meters if the RasterLayer is not projected (+proj=longlat) and in map units (typically also meters) when it is projected.
popcenterdist=distance(popcenter)
plot(popcenterdist)
Will sea level rise affect any major population centers?
Steps:
popcenter to resolution of dem using method=ngbpopcenter areas that flood according to flood2.Will sea level rise affect any major population centers?
popcenter2=resample(popcenter,dem,method="ngb")
floodpop2= flood2==1 & popcenter2
floodpop2=mask(floodpop2,floodpop2,maskval=0)
plot(flood2);plot(floodpop2,add=T,col="red",legend=F);plot(bgd,add=T)
vpop=rasterToPolygons(popcenter, dissolve=TRUE)
gplot(dem,max=1e5)+geom_tile(aes(fill=value))+
scale_fill_gradientn(
colours=c("red","yellow","grey30","grey20","grey10"),
trans="log1p",breaks= log_breaks(n = 5, base = 10)(c(1, 1e3)))+
coord_equal(ylim=c(21,25))+
geom_path(data=fortify(bgd),aes(x=long,y=lat,order=order,group=group),size=.5)+
geom_path(data=fortify(vpop),aes(x=long,y=lat,order=order,group=group),size=1,col="green")
Warning: very slow on large rasters…
Uses rgl library.
plot3D(dem)
decorate3d()
50 different styles illustrated here.
Overlay population with drape
plot3D(dem,drape=pop3, zfac=1)
decorate3d()
distance() can be slow!